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04 ^ The center-of-mass dynamics of a vortex in the surface region of a Bose-Einstein condensate is 

investigated both analytically using a variational calculation and numerically by solving the time- 

^ ' dependent Gross-Pitaevskii equation. We find, in agreement with previous works, that away from 

the Thomas- Fermi surface, the vortex moves parallel to the surface of the condensate with a constant 

velocity. We obtain an expression for this velocity in terms of the distance of the vortex core from 

("^ ' the Thomas- Fermi surface that fits accurately with the numerical results. We find also that, coupled 

^Sl ' to its motion parallel to the surface, the vortex oscillates along the direction normal to the surface 

around a minimum point of an effective potential. 



i-C ' PACS numbers: 03.75.Fi, 67.40.-w, 32.80.Pi 
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I. INTRODUCTION 

-(— > , 

^_> I Following the experimental observation of vortices in Bose-Einstein condensates m , [ a [j jj^, M, 0, S, Q , vortex 

. • dynamics has been a subject of intensive interest [13, [O, ]12, Ej Hj llM llfl ^2, lH) Um 120, UM, ^^ ■ Theoretical studies 

"trt ■ have been performed on the center-of-mass motion of the vortex core in the bulk of the condensate [lO, \lM, l23 |. 

a References flOJ, I24I have studied the dynamics of the core taking into account the dissipation from the thermal 

cloud. Dissipation causes the vortex to spiral out of the condensate giving it a finite lifetime in agreement with 
'^ , the experimental observations [2]- Once they reach the surface of the condensate, the vortices decay into surface 
Ph ■ excitations. It is also believed that vortices enter the condensate from the surface at the point when surface excitations 
O " become unstable [10, ^M, 123, 0, 123, 120, 123 • Hence, it is worth investigating the details of the dynamics of a vortex 
in the surface region of the condensate. 

E. Lund et al. [2a | have pointed out that near the surface of the condensate, the quadratic trapping potential can 
O^ ■ be approximated by a potential that is linear in distances normal to the surface. As a result of this approximation, 
^ ' the equilibrium density profile will be also linear in distances normal to the surface of the condensate except for 
distances very close to the Thomas-Fermi surface. An immediate question then arises: How will the vortex move in 
such an inhomogeneous density background? The fact that the density is linear makes it possible to perform analytical 
vQ treatment to the problem provided that the size of the core of the vortex is much less than the distance over which the 
f~^ ' background density varies significantly. This question has been addressed extensively by Anglin [20] , where the author 
■TJ" [ used a boundary-layer theory to calculate both the density profile and phase of the vortex, from which he obtained 
^D ■ the equilibrium velocity of the vortex. It turned out that the component of the velocity in the direction normal to 
the surface of the condensate vanishes, while the component parallel to the surface is a constant that depends on the 
distance between the core of the vortex and the surface. In this paper we address basically the same problem but using 
a variational calculation. The main aim is to demonstrate that the variational calculation performed here generates 
the previous knowledge about the vortex motion near the surface and captures the main features of it. We go one 
step ahead of previous works by obtaining a rather more accurate expression for the component of the vortex velocity 
Q ' parallel to the surface. Furthermore, we show that, within few coherence lengths from the Thomas-Fermi surface, the 
CJ I vortex oscillates along the direction normal to the surface around a minimum point of an effective potential. These 
oscillations are coupled to its motion parallel to the surface and disappear when the vortex is away from the surface. 
In our variational calculation we use an ansatz wavefunction that takes into account the equilibrium density profile 
and phase of both the background and the vortex. The variational parameters are the coordinates of the vortex 



H ] core and the center-of-mass velocity components normal to the surface of the condensate and parallel to it. The 
equations of motion for these parameters are then derived from a lagrangian that corresponds to the time-dependent 
Gross-Pitaevskii equation. 

Using a linear density profile for the background is an approximation that breaks down at the Thomas-Fermi surface. 
Our variational calculation is, thus, expected to be accurate only away from the Thomas-Fermi surface. Furthermore, 
we neglect the internal dynamics of the core of the vortex. This is also justified away from the Thomas-Fermi surface 
where the size of the core of the vortex is almost constant. 

In the following, we start, in subsection III Al by writing down the ansatz wavefunction. Using this wavefunction, we 
calculate the lagrangian from which we derive the equilibrium properties in subsection lll Bl and the equations of motion 
for the variational parameters in subsection III CI In Section IIIII we present the results of the numerical simulation 



of the vortex motion, and compare them with those of the variational calculation. In Section IIVI we summarize our 
main conclusions. 



II. VARIATIONAL APPROACH 

A. The Lagrangian 

The dynamics of the time-dependent Bose-Einstein condensate is described by the Gross-Pitaevskii equation 
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Here g is the effective two-particle interaction which is proportional to the s-wave scattering length a according to 
g = AnaTi /m, where m is the mass of an atom, and /i is the equilibrium chemical potential. The harmonic trapping 
potential ^(r') is given by 



V{v') = -mujnr 
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where ujq is the characteristic frequency of the trap. Here, we take for simplicity V{r') to be isotropic. For distances 
near the surface of the condensate, the quadratic trapping potential V(r') can be approximated by a linear potential 

MM 



V{r') ~ V{R) + Fx , 



(3) 



where a; is a coordinate normal to the surface of the condensate, R is the Thomas-Fermi radius given by V{R) = fi, 
and F — niujgR is a force constant. The last equation is obtained by using the transformation r' = r -I- Rx in Eq. ([2]). 
This suggests a planner geometry in which the surface of the condensate is approximated by an infinite plane. The 
x-coordinate will be normal to the plane while the y and z coordinates will be parallel to it. The origin of this 
coordinate system is then shifted from the center of the condensate to the surface. As a result of this shift, the bulk 
of the condensate occupies the region a; < and the surface region will be near x — 0. This is illustrated in Fig. [H 
In this planner geometry, the Gross-Pitaevskii equation takes the form 
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The equilibrium density profile riQix) — \ip{x)\'^ is obtained from the last equation by setting the time derivative to 
zero 
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For large negative values of x, the kinetic energy operator can be neglected |28l . |30| leading to the Thomas-Fermi 
equilibrium density nj^(a:) — —{F/g)x. 

We consider a vortex with a core located at r = xo{t)x + ya{t)y and axis along the z-direction as shown in Fig. [TJ 
The vortex has a center-of-mass velocity v — Vx{t)Si. + Vy(t)y. Our trial wavefunction is given by 
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X{p) exp («$(p, (j))) exp [i{m/K){xVx + yvy)] . 



(6) 



Here, the first factor on the right hand side is y^z^^a?). The second factor corresponds to the square root of the 
density of a vortex in a uniform background. The function $(p, 0) in the exponential is the phase of the vortex 
that takes into account the inhomogeneouty of the background density. The last factor represents the center-of-mass 
velocity of the vortex. The polar coordinates p — ^J{x — XQ{t)Y + {v — yoit))"^ and (j) = tan~^ {{x — xq) / {y — yo)) 
have the core of the vortex as their origin, as indicated in Fig. [1] Near the surface, the Thomas-Fermi approximation 
breaks down and the density background is nonlinear. The kinetic energy starts to become important leading to an 
exponential decay of the density. Our trial wavefunction does not take this nonlinearity into account. Therefore, the 
results of this variational calculation are not expected to be accurate in the surface region where the Thomas-Fermi 
approximation breaks down. 



For the phase of the vortex we use an expression similar to the one derived by Anghn [2C 
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The first term of this expression represents the usual (/)-dependance of a vortex in a uniform background. The second 
term results from the linear inhomogeneouty of the background. It should be noted that the angle cf) defined in this 
paper is the complement of that defined in Ref . [20] . This expression is obtained using a boundary-layer theory which 
is described here briefly. Vo rtices have structure in their density over a small length scale of the order of the coherence 
length ^0 = l/-\/87rano(a;o), and phase that extends over the whole system. Here, ^o is the coherence length calculated 
with the background density at xq, and S = {h/mFY'^ is the surface depth [2»|. In the boundary-layer theory, the 
vortex problem is first solved exactly in the mner region p ^ £,o and the owier region p S> ^o- The inner asymptotics of 
the outer solution and the outer asymptotics of the inner solution are then obtained by expanding the outer and the 
inner solutions in powers of two perturbation parameters characteristic to the two regions. Finally, the two solutions 
are matched to determine the solution in the intermediate region. The velocity and phase of the vortex are obtained 
as a result of this matching procedure. The phase, given by Eq. ([7]), corresponds to the inner asymptote of the outer 
solution. Thus, this expression is valid only near the vortex and is not expected to be accurate for p ^ ^o- However, 
in the present calculation, we take this expression to be valid for all values of p. We account for the fact that it is 
not valid for distances away from the vortex core, by introducing in the argument of the logarithm the parameter a, 
which works as an effective cutoff on the logarithmic term. The value of a is then determined by fitting the velocity 
calculated variationally in subsection IIIBI with the one obtained from the numerical solution of the time-dependent 
Gross-Pitaevskii equation in section Hill 

The vortex density profile is taken as the square of [3l| 
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It turns out in subsection III CI that the equations of motion of yo(i) and Vy(t) are independent of x- Therefore, 
specifying the functional form of x is not necessary for the y-component of the motion. 

Due to the nonlocal nature of the vortex excitation, the energy of the vortex diverges logarithmically with the size 
of the system [S^l- To be able to perform the variational calculation for a finite and constant number of atoms N, we 
consider in calculating the lagrangian only the atoms within a cylinder of radius b and length I such that the cylinder 
axis coincides with the vortex axis. The fact that we restrict the number of atoms within the cylinder to be constant, 
is guaranteed by requiring ip{r,t) to be normalized to N. The normalized wavefunction takes the form 
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where 



Ni= / dpp / dcf>x^- 
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The lagrangian that corresponds to Eq. ^ is given by 
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where AE[ip, Tp*] is the energy associated with the presence of the vortex. This energy is obtained by subtracting the 
energy of the vortex-free background from the vortex energy |32l | , namely 
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and the vortex-free energy is 
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The terms on the right hand side of the last equation represent the mean-field and trapping potential energies of the 
vortex-free background, respectively. These were calculated using the wavefunction 
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The prefactor of -^vf guarantees its normalization to N . For the subtraction of the background energy from the vortex 
energy to be correct, it is essential that the number of atoms of both systems to be the same. 

Using the trial wavefunction Eq. ^ in the lagrangian of Eq. (fTTj) . we show in Appendix \K\ that the lagrangian per 
atom takes the form 
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where the coefficients A^i, .. .,N^ are functions of b, and N2 and A^4 have also Xq dependence ( See Appendix [XI). 
The dot on j/Qj ^x, and Vy denotes derivative with respect to time. The parameter 7 — gN/Nil is an energy that 
characterizes the average mean- field energy since N /Nil represents the average density within the cylinderical volume 
surrounding the vortex. This lagrangian will be the basis of the calculations of the equilibrium and nonequilibrium 
properties in the next two subsections. 



B. Equilibrium Properties 



The equilibrium values of the variational parameters can be obtained by minimizing the energy functional with 
respect to the variational parameters. The energy functional can be readily obtained from Eq. (J16p by setting time 
derivatives to zero and then multiplying by -1, namely 
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Minimizing A_B with respect to Vx gives Vx — Q. As we shall see in the next subsection, Vx is proportional to xq in 
the limit xq ^ —5. Therefore, Xq will be a constant throughout the motion of the vortex. In other words, the vortex 
will be moving parallel to the surface of the condensate in agreement with previous works [la . |20| and the results of 
the numerical simulation of section IIIII Minimizing the energy with respect to Vy yields 
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In the hydrodynamic limit, 6 3> C, the coefficient A'4(xo) is expanded in powers of ^/b as A'4 
O {{£,/b)^ In (6/C))- Using this expression in the last equation, it takes the form 
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The quantity h/2mxQ is the equilibrium velocity of a vortex in a homogeneous density background located at a distance 
xq from a hard wall and is moving parallel to it. The logarithmic term corresponds, therefore, to the contribution of 



the linear inhomogeneouty in the background density. The value of ab is determined by fitting this expression to the 
numerical values of Vy{t). This is performed in section Hill where it turns out that ab f^ 0.11. Thus Eq. P^ reads 
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This is to be compared with the expression derived by Anglin [2 
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using a boundary-layer theory, namely Vy{t) = 



{Ti/2mxo{t)) (1.96 + (3/2) In {xo/6)). In Fig.[2l we plot Vy as a function of a;o using Eq. ((20)) . the expression of Anglin, 
and the numerical solution of the Gross-Pitaevskii equation. This figure shows that while the expression of Anglin is 
accurate for large values of xq, it deviates from the numerical values of Vy for values of xq of order 6. On the other 
hand, Eq. ((20)) describes accurately Vy for both large and small values oi xq. Equation ([20]) represents one of the main 
results of this paper as it shows that the variational calculation does indeed lead to a rather accurate description of 
the vortex velocity parallel to the surface of the condensate. 

Substituting the equilibrium expression of Vy from Eq. (jlSp in Eq. (|17p . we obtain an energy functional that is a 
function of Vx and xq , namely 



AE[vx,xo]/N = -mvl + U{xo) 
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where 
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is a function of xq that can be considered as an effective potential for the vortex. Here, we have neglected a term 
containing Nf since it is one order of magnitude smaller than the logarithmic term of Eq. ^ . Inspection shows that 
this effective potential has a minimum for negative values of Xq . This means that the vortex can acquire oscillations 
around this minimum. Taking this into consideration and the fact that for large xq , we have Vx = io, suggests following 
a collective coordinates approach to find the frequency of oscillation [29|. In this approach the energy functional, 
Eq. ([2T|) . is put in the form of that of a simple harmonic oscillator /^E[vx,xq\ = mv'^/2 + k{xQ — Xeq)^/2, where Xoq is 
the value of xq at which U{xif) has a minimum. The frequency of oscillation will be simply to = yjkjm. This is achieved 
by expanding the coefficients iVi, N2{xq)^ N^, Nq, and A^7 in powers oi£,/b keeping terms up to the second order, namely: 
N2{xo) = (1 + ^l-{b/xoy)/2b^ + O ((e/teo)2) , N^ = by 2 + 2e In (&/0 + H\Uhf In^ (fe/O - 
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Expanding U{xo) around Xeq, we obtain 
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The energy functional becomes 



AE[xo,Vx]/N = T^ln- 

The frequency of vortex oscillation is then given by 
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Noting that j — h /2?Tif^ ~ n(xo) and that in the limit |a:o| 3> S, the Thomas-Fermi approximation gives n{xo) ~ |a:o|, 
we conclude that tu ^ y/\xf)\~'^lnb/(,. Thus, for large Xq, these vortex oscillations disappear and Xq will be constant 
as we have also found above by minimizing AE with respect to Vx- 

In the next subsection, we derive equations of motion that describe the vortex motion and oscillations. 



C. Equations of Motion 



The equations of motion follow from the Euler-Lagrange equations of the lagrangian Eq. (|16p. namely 
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The prime on N2{xq) denotes a derivative with respect to a;o. In obtaining the last equation in terms of a;o only, we 
have used Eqs. ([29|) and (f30|) . In addition, we have neglected terms containing Nl{xo) and A'4(a;o)A^4(a;o) since these 
terms are of an order of magnitude less than the logarithmic term in the phase <&(/o, (j)) of Eq. ([7]). 
From the first two equations, we obtain quite generally 



i/o(0 = o. 
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This is a general result in the sense that it does not depend on the functional form of x, the upper limit of the radial 
integration b, the cutoff of the logarithmic part of the phase a, or xq- The coordinates of the vortex are then described 
completely by Eqs. dSJ) and ((5^ . 

In the limit |a;o| ^ S, Eqs. ([50)1 and ([?!]) reduce to v^ = xq and xq — 0, respectively. Therefore, xq will be in 
general a linear function of i, namely xo{t) — cit + C2, where ci and C2 are constants. We have shown in the previous 
subsection that Vx — for the energy to be minimum. This results in that ci — 0, and thus a;o is a constant as we have 
also found in the previous subsection. We notice also that, in this limit, Eq. ([29|l reduces to Vy ^ yo. This implies 
that in the limit |xo| 3> S, the variational parameters Vx and Vy are the conjugate velocities of xq and j/Qj respectively. 

In the limit 6 3> C, the coefficients iVi, . . . ^N-j can be expanded in powers of ^/6, as was shown in the previous 
subsection, and Eq. ((3T|) simplifies to 
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where U{xo) is the effective potential of the vortex defined in Eq. ((23|) . Expanding the last equation around 
using Eq. (^5]) . we get 
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Using the initial conditions xo(0) — Xoq + Xq and io(0) = 0, the solution of the last equation reads 
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where fl = ^/{^ — &^/4x^„) and w is given by Eq. ((27|) . Notice that in the limit b ^ 6, the quantity b'^/x'^„ — > 0, 
which leads to O — > w. Equation (|3ip can be solved numerically for given values of b. We plot this numerical solution 
in Fig. [Altogether with analytical solution Eq. ([55]) . This figure shows that the analytical solution is a reasonable 
approximation for the exact solution of Eq. ([5T|) . 

To ensure that these oscillations are not an artifact of our variational calculation, we have solved numerically the 
time-dependent Gross-Pitaevskii equation for initial vortex distances from the Thomas-Fermis surface that are of 
order d. We found that the vortex indeed oscillates in the x-direction as predicted by the variational calculation. 
Furthermore, we show next that we can obtain good agreement between the variational calculation and the numerical 



one for both the frequency and amphtude of the oscillations. In the present variational calculation, the amplitude 
of vortex oscillation X^ and frequency uj are essentially determined by Eqs. ([M]) and ([?7|) . Scaling, length to 6 and 
energy to h jvnS^ = F5, these two equations reduce, in the zeroth order of ^/6, to: x^q — {S/S,)'^S and uj/{F/m6) — 
i^/S)'^ In {b/^), where we have used fi /2m^'^ — 7. Thus, to obtain the amplitude in units of 5 and the frequency in 
units of Fjmb, we need only to specify ^/(5 and hjb. The parameter hjb is a free fitting parameter that will be set 
a value to get the best agreement with the numerical solution. In Fig. [31 we plot a;o(t) obtained from the numerical 
solution of the time-dependent Gross-Pitaevskii equation together with the solution of Eq. (j3ip using the same initial 
conditions. The vortex is initially located at Xg = —3.1(5 from the Thomas-Fermi surface and then evolved starting 
from rest. The best agreement between the variational solution and the numerical one was obtained for h ~ hb. The 
coherence length associated with this figure is calculated by fitting the numerical density profile to Eq. ([5]), using 
^ as a fitting parameter. This was done for both the x- and y-cross sections of the vortex density profile as shown 
in Fig. m The average value of ^ turns out to be ^ ~ 0.425. This is consistent with the Thomas-Fermi estimate 
^ = 2m'y /% ~ 0.4(5, were we have used 7 = gn(xo) = ^Fxo- Since the coherence length obtained using the Thomas- 
Fermi approximation is approximately equal to the one obtained numerically, we conclude that during the whole 
time interval of Figs. [3] and IH the vortex is moving in a region where the Thomas-Fermi approximation is valid. At 
first instance, we may take the average value of the coherence length, namely ^ = 0A6, and then substitute it in 
Xcq and Lo. This will lead to Xcq — —6(5. However, Fig. [3] indicates that the vortex oscillates around Xcq — —2.8(5. 
This discrepancy originates probably from the difference in geometry of the system in the variational calculation and 
the one in numerical calculation. In the variational calculation, the vortex is moving in an infinite background, but 
the energy integrations are performed within a cylinder of radius b — 5(5. Furthermore, the number of atoms within 
the cylinder is constant. In the numerical calculation, the system is a rectangular grid of dimensions 86 x 40(5. The 
number of atoms within the whole grid is constant. Therefore, while in the variational calculation the number of 
atoms is kept constant within an area of 7r&^ ~ 80(5^, in the numerical calculation the number is kept constant within 
an area of 320(5^. Thus, there is a difference in the density which leads to a difference in the coherence length. In 
other words, the coherence length values shown in Fig. [1] are not exactly those which we should use in the variational 
calculation to seek agreement between the results of the two calculations. We found, as shown in Fig. [31 that the best 
agreement with the numerical values is obtained using ^ = 0.86. 

We show also in Fig [31 vortex trajectories for initial distances xo(0) = —4(5 and a;o(0) — —5(5. There is a transient 
behavior in the beginning before acquiring the oscillatory behavior. For initial starting distances larger than 7(5, the 
oscillations disappear and the vortex moves in a straight path parallel to the surface as predicted by the variational 
calculation. 

III. NUMERICAL SOLUTION 

In this section, we describe our numerical procedure in solving the time-dependent Gross-Pitaevskii equation. 
Scaling length to 6, energy to h /m6'^, and time to m6^/h, Eq. ^ takes the following dimensionless form 
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Mr,t) = i—Mr^i) , (36) 
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where the tilde denotes scaled quantities and |'0P is scaled to l/47ra(5^. It is interesting to notice here that the 
scattering length disappears as a result of the scaling we employ here. 
We start by solving the time-independent Gross-Pitaevskii equation 
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with a vortex-free initial state. We use the imaginary-time evolution method to solve this equation. The resulting 
density no{x) is the background density profile for the vortex. We multiply this background density by a density 
profile of a vortex in a uniform background, namely x^ given by Eq. ([8]) . The resulting density profile can not yet be 
used as an initial state for solving the time-dependent Gross-Pitaevskii equation. This is so since x of Eq. ^ is not the 
exact density profile of the vortex. When we simulated the dynamics of the vortex using this density profile, the core 
of the vortex acquired breathing oscillations during the vortex motion. These oscillations produced circular density 
waves emitted out of the vortex core through the background and then refiected from the boundaries to interfere with 
vortex affecting its trajectory. Since we are interested only in the center-of-mass motion of the vortex, it is necessary 
to use the exact density profile of the vortex to avoid exciting such density waves. This can be obtained by evolving 
in imaginary time for a certain time period the density nox^- In other words, we resolve Eq. (|37p with -^/noX^ as 
an initial wavefunction. The resulting profile contains the exact background and vortex core structure. This profile 
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is then used as the initial profile for the time-dependent Gross-Pitaevskii equation. We found that density waves 
emerging from the vortex core have now almost disappeared. However, evolving in imaginary time results in the 
vortex sliding towards the surface of the condensate and ultimately disappearing there. In the present procedure, we 
evolve for only a finite time interval long enough to ensure adequate suppress of the density waves before the vortex 
reaches the surface. The distance from the surface reached by the vortex at the end of this time interval is to be 
considered as the initial distance for the real time evolution. We checked that our results for the real-time dynamics 
of the vortex are not affected by the length of this time interval once the density waves are significantly suppressed. 

We employed a fourth-order Runge-Kutta method to integrate the time derivative of the time-dependent Gross- 
Pitaevskii equation. The spacial grids we used are 86 x 40(5 and 85 x 80i5 with spacial descretization of Ax — 0.16 
and time descretization At = 10~^. The long side of the grids is along the y-axis since the vortex is faster in that 
direction. We have used the boundary conditions that the derivative of density with respect to y vanishes at the two 
edges of the grid normal to the y-axis. The other two edges were left free. We have used larger grids to check that 
the boundaries do not affect the results obtained here and found that the above-mentioned grid sizes are adequate for 
that purpose. We have also calculated the total number of atoms and energy to find that they are fluctuating within 
less than 0.1%. The fluctuations were random and not correlated to the vortex motion. 

To test the accuracy of our numerical calculation, we start by simulating the dynamics of a vortex near a hard wall 
in a condensate of uniform background density. The vortex velocity is calculated theoretically to be Vy = h/2mxQ 
[2(j|. In Fig. [51 we show that our code accurately generates such a dependence for the vortex velocity. 

Having checked the accuracy of our numerical procedure, we simulate the dynamics of the vortex in the surface 
region of the condensate. In Fig. [51 we show the results in the form of a number of snapshots of the density of the 
condensate visualizing the realtime motion of the vortex core along the surface of the condensate. In Figs. [3] and 
m we plot the vortex coordinates Xq and j/o versus time. These two figures show that while the vortex is moving 
with a constant velocity in the y-direction, it oscillates in the x-direction as predicted in the previous section. To 
investigate the xo-dependance of Vy, we calculate the numerical value of Vy for a number of different values of xq. The 
results are shown in Fig. [2] where we also plot the theoretical prediction of Anglin [2d| and the result of our variational 
calculation. The formula derived by Anglin is supposed to be valid away from the Thomas-Fermi surface. Figure [2] 
shows indeed that while this formula agrees well with the numerical calculation for large xq, it deviates from it for 
values of xq of order 6. On the other hand, this figure shows that our formula, derived in subsection IIIBl is in a 
better agreement with the numerical data especially for xq < 106 where the formula of Anglin starts to deviate from 
the numerical values. 

We have also calculated the coherence length along both the x and y-directions, which is of the order of the size of 
the core of the vortex. This is plotted in Fig. [H This figure shows that the core sizes change slightly over the whole 
time interval considered. The average value of ^ in this figure agrees with the Thomas- Fermi estimate as explained 
in subsection III CI In the previous section, we performed the calculation using the Thomas-Fermi approximation for 
the background density and neglected the dynamics of the size of the core. It is clear now that this is indeed justified 
here. 

IV. CONCLUSIONS 

We have shown that our variational calculation describes at least qualitatively the motion of a vortex near the 
surface of a Bose-Einstein condensate, and generates results that agree with previous calculations [l^ \1^. These 
qualitative results were then supported and quantified using a numerical simulation of the vortex motion. 

One of the main findings of this work is Eq. ([20)l which gives the y-componcnt (parallel to the surface) of the 
vortex velocity as a function of it is distance from the Thomas- Fermi surface xq. It agrees qualitatively with that 
of Anglin f20J , and describes more accurately Vy for values of Xq of order of 6 as shown in Fig. [21 This figure shows 
clearly that Eq. ([50)1 can be used for both small and large xo/6, unlike the result of Anglin |20| which is accurate only 
for large xo/6. 

Another main result of this work is that, if the vortex is moving within few coherence lengths ^ from the Thomas- 
Fermi surface, it will have an oscillatory component of its velocity in the direction normal to the surface. These 
oscillations disappear for xq <S^ —6 and the vortex will have only the y-component of the velocity. To further understand 
the physics of these oscillations, we have shown in subsection IIIBl that the vortex feels an effective potential that 
is a function of xq and has local minimum for negative values of xq . This effective potential is composed mainly of 
a sum of a mean-field-energy and a trapping potential energy part as shown in Eq. (|22p. These two parts compete 
such that the effective potential has a minimum at Xoq — —((5/^)^(5. Using a collective coordinates approach, we have 
shown that the vortex oscillates around this minimum with a frequency u; = y''{F'^£^'^ /8mj^) In (b/^). As discussed 
in subsection IIIBl this frequency is of order less than -y/lxol"^ In (|a;o|/^). Both this frequency and the xo-dependent 
part of U{xo) vanish for large xq. This means that the vortex oscillations will disappear away from the surface of the 



condensate and will be noticeable only when ^/xq is not negligible. To the best of our knowledge, these oscillations are 
predicted here for the first time. As a possible alternative explanation to these vortex oscillations near the surface, we 
mention the possibility that the vortex being dragged by a collective or surface wave. In our numerical procedure, the 
evolution in imaginary time which generates the initial state for the real time evolution, may include some excitations 
of collective modes. Further investigation is needed to verify this argument. 

It was shown in Ref. [20] that the energy of a vortex near the Thomas-Fermi surface is not infrared divergent as 
in the bulk. This is so since the vortex is localized to within its distance from the surface. To avoid these infrared 
divergencies in the energy, we have introduced an infrared cut-off on the energy. This is equivalent to using the 
infrared-free flow pattern of Ref. [20| which gives finite energy. 



APPENDIX A: DETAILS OF THE CALCULATION OF THE LAGRANGIAN 



In this Appendix we present the details of calculating the lagrangian Eq. (|16p . This is performed by substituting 
the trial trial wave function Eq. ^ in Eq. pTjl . 

We start by calculating the energy functional Eq. (fT3)) . According to the geometry illustrated in Fig.[Tl the cartesian 
coordinates x' = x — xq and y' = y~yo with origin at the core of the vortex are related to the polar coordinates p{x' , y') 
and 4>{x' ,y') by p{x',y') — \/ x'"^ + y'"^ and 4){x' ,y') = ia.n~^{y' /x'). This leads to that d(t>/dx — d(l)/dx' — cos(/)/p 
and d(j)/dy = d(j)/dy' = —sincp/p. In addition, we have dp/dx = dp/dx' = smcf) and dp/dy — dp/dy' — coscj). The 
derivative of the trial wave function ip{x, y, t) with respect to x, can thus be shown to take the form 



dx 



The contribution to the kinetic energy is obtained by integrating \dil] /dx]"^ with respect to p and i 
reduces to 
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where Xp denotes dx{p)/dp. It should be noted that we neglected here the contributions of order 

that originate from the logarithmic term of the phase $(p, 0). In deriving this term, Anglin [2C 

turbative expansion of the phase of the vortex near the core using ^/xq as the perturbation parameter, where S, is 

the coherence length defined in subsection III Al Since in Anglin's calculation only terms of order {£,/ p) In (ap/xo) are 

kept, it will be inconsistent if we keep here higher order terms. We keep, however, the higher order terms that do not 

originate from the logarithmic term of the phase. 

Similarly, for the derivative with respect to y, we have 
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which leads to 



2ir 



The kinetic energy is then given by 
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The potential energy associated with the external force is 
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The mean-field energy is 
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Combining these three contributions, we find the energy functional per atom 
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where 7 = gN/Nil and we have defined the coefficients 
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The relevant energy is the energy change due to the presence of the vortex. Therefore we have to subtract from the 
above energy functional the energy of a vortex-free background 
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Thus, the energy difference associated with vortex reads 
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The time derivative of the trial wave function is given by 
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(All) 
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Using that dp/dt — [dp/dx'){dx' /dt) + {dp/dy'){dy' /dt) = — ipsinc/) — y'^coscj) and similarly d(j)/dt — {—x'qCOScJ) + 
ijQ sin (j))/p and Eq. ([7]), the part of the lagrangian (per atom) containing the time derivatives takes the form 

- nf^il - 2Ni{xQ)) - m{xov, + yovy) - Ns—^ . (A12) 

Axq ZXq 

The lagrangian per atom, Eq. (J16p . is then obtained. 
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FIG. 1: The planner geometry and coordinate system used in this paper. The smaU circle represents the core of the vortex 
which is of the order of the coherence length ^. The coordinates of the core of the vortex are xq and j/o- The radius of the 
condensate is 7?. The bulk of the condensate is located in the region a; < 0. The point P has cartesian coordinates x' and y' 
related to the polar coordinates p and by x' = psini;/) and y' = pcostf). 
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FIG. 2: The vortex velocity in the y-direction as a function of the magnitude of the distance between the vortex core and the 
Thomas-Fermi surface. The velocity is in units of 5/t. The surface depth 5 is defined by h^ /mS'^ — F5, and the time r is 
defined by r = h/FS. The dots on the upper curve are the results of the numerical solution of a vortex moving near the surface 
of the condensate. The dots on the lower curve correspond to a vortex moving parallel to a hard wall in a uniform density 
background. The solid curve is the result of the variational calculation of section [III given by Eq. (|20|l . The long-dashed curve 
corresponds to the expression derived by Anglin [2^]. The short-dashed curve is given by Vy — h/2mxQ. 
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FIG. 3: The x-coordinate of the vortex core measured from the Thomas-Fcrmi surface. The points are the results of the 
numerical solution of the time-dependent Gross-Pitaevskii equation, Eq. ^ for different initial distances. The initial velocity 
is taken to be zero. The solid and dashed curves are the results of the variational calculation of subsection III CI The solid 
curve is the solution of Eq. pi|l . and the dashed curve is an analytic approximate solution given by Eq. p5[) . The value of b 
used here is 55. 
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FIG. 4: The coherence length associated with the x- and j/-cross sections of the vortex density profile. It is calculated by fitting 
the numerical density profile of the vortex to Eq. (|8j . Circles represent the coherence length associated -with the y-cross section 
of the vortex density profile, and diamonds represent the coherence length for the a;-cross section. 
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FIG. 5: Density profile of tfie vortex. Tfie upper row of pictures corresponds to tfie initial state with the vortex core located 
at a;o(0) = —3.15. The left column of pictures is a front view in which the x-axis is normal to the picture towards the viewer. 
The right column is a side view in which the j/-axis is normal to the picture towards the viewer. The vertical line in the right 
column of pictures is to reference the core location with respect to the initial position, and shows clearly that the vortex core 
oscillates in the a;-direction. 
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FIG. 6: The vortex j/-coordinate. Starting from the upper curve, the values of the initial position of the vortex core xq are: 
—35, —4(5, -5(5, and —75. The constant slope indicates that the velocity component parallel to the surface is constant. 



